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Abstract.  The  direct  simulation  Monte  Carlo  (DSMC)  method  is  used  in  a  time- dependent  manner  to  simulate  three- 
dimensional  micro  Couette  flows.  An  artificial  forcing  that  mimics  the  centrifugal  force  in  the  Taylor  problem  has  been 
applied  to  the  flow.  The  sampled  behaviors  of  the  resulting  flow,  including  the  averaged  properties  and  disturbances,  are 
studied.  The  computations  have  been  performed  using  a  parallel  computer  cluster.  The  results  presented  include  those 
with  various  channel  heights,  plate  speeds,  and  the  forcing  level.  These  changes  also  result  in  changes  in  the  flow 
Reynolds  number  and  Knudsen  number.  Spatially  coherent  flow  patterns  can  be  identified  in  the  averaged  flow  and  the 
disturbance  flow  fields.  The  results  indicate  that  the  discrete  approach  can  capture  unsteady,  three-dimensional  vortical 
flow  structures.  In  cases  with  strong  forcing,  the  disturbance  energy  spectra  show  significant  content  above  statistical 
scatter. 


INTRODUCTION 

Microelectromechanical  systems  (MEMS)  technologies  have  had  significant  impacts  in  different  areas  such  as 
biosciences,  computer  sciences,  and  telecommunication.  Many  of  these  realized  devices  use  fluid  as  a  working 
media  and  their  designs  have  often  been  guided  by  correlations  derived  from  the  fluid  flow  behavior  observed  at 
macroscales.  With  the  increasing  demand  for  higher  system  complexity  and  performance,  it  is  important  to  develop 
a  better  understanding  of  the  fluid  flow  characteristics  at  the  microscale. 

Many  recent  experiments  in  fluidic  microchannels  have  reported  significant  differences  in  the  heat  and 
momentum  transfer  coefficients  compared  with  those  at  the  macroscale.  For  example,  while  the  friction  factor 
varies  inversely  with  the  Reynolds  number  when  the  Reynolds  number  is  small,  the  proportionality  constant  does 
not  agree  with  the  conventional  correlation.  The  friction  factors  are  also  found  to  diverge  from  this  inverse 
proportionality,  which  at  macroscale  indicates  a  change  of  flow  characteristics  from  that  of  a  laminar  flow  to  a 
turbulent  flow,  at  smaller  Reynolds  numbers  than  those  commonly  observed  in  the  corresponding  large  channels. 
This  early  change  of  flow  characteristics  has  been  attributed  to  the  effects  of,  for  instances,  gas  rarefaction  and  other 
surface  mechanisms,  the  large  change  of  the  flow  Reynolds  number  along  the  microchannels,  and  the  likely 
experimental  uncertainties  in  microscale  measurements.  Advanced  measurement  techniques  such  as  molecular 
tagging  can  produce  more  detailed  quantitative  data  for  microflows,  which  may  eventually  lead  to  a  better 
understanding  of  the  apparent  microflow  transition.  In  this  paper,  results  obtained  in  recent  micro  Couette  gas  flow 
simulations  will  be  presented.  The  DSMC  method  pioneered  by  Bird  [1]  is  used.  With  an  artificially  applied, 
adjustable  forcing,  the  flow  model  is  used  to  study  capturing  the  unsteady,  three-dimensional  flow  disturbances  in 
microflows. 

The  DSMC  method  has  been  applied  to  many  different  fluid  dynamics  problems.  It  has  been  widely  used  in  the 
simulations  of  rarefied  gas  flow.  The  method  has  also  recently  been  applied  to  calculate  the  fluid  flow  and  heat 
transfer  behaviors  of  microchannels. [2,3]  DSMC  has  been  employed  to  study  the  low-density  limits  of  a  number  of 
flow  instabilities.  The  centrifugal  instabilities  in  Taylor-Couette  flows  have  received  the  most  attention. [4,5]  The 
formation  of  Taylor  vortices  were  clearly  demonstrated  for  a  range  of  Knudsen  number  in  these  two-dimensional 
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(axisymmetric)  flow  computations.  The  convective  instability  associated  with  the  Rayleigh-Benard  cells  has  also 
been  studied.  [6] 

DSMC  applies  a  molecular  description  of  gases.  Gases  are  recognized  as  a  myriad  of  discrete  molecules  and  the 
positions,  velocities,  and  states  of  molecules  are  recorded  at  all  time.  The  computational  molecules  are  allowed  to 
move  and  collide  in  the  physical  space.  Macroscale  flow  properties  can  be  obtained  by  sample  averaging.  For  gas 
flows  in  microscaled  geometry,  the  Knudsen  number  based  on  the  characteristic  size  of  the  geometry  can  be  large. 
The  applicability  of  the  linear  constitutive  relations  in  the  Navier-Stokes  equations,  which  are  used  in  the  direct 
numerical  simulations  of  macroscale  flows,  becomes  questionable  due  to  the  large  flow  gradients.  In  fact,  when  one 
considers  the  local,  smaller  scale  flow  events,  the  value  of  the  local  Knudsen  number  can  further  increase.  The 
DSMC  method  is  based  on  first  principles  and  is  not  limited  to  small  Knudsen  number.  DSMC  is  appropriate  for 
the  direct  numerical  simulations  of  such  microflows.  As  the  DSMC  method  uses  a  molecular  description  of  gases, 
there  have  been  concerns [7]  whether  DSMC  is  capable  of  resolving  vortical  fluid  motions  due  to  the  lack  of  an  exact 
conservation  of  angular  momentum  in  collisions.  Results  obtained  by  Nanbu  et  al.[8]  shows  that  as  long  as  the  cell 
size  is  sufficiently  small,  the  total  angular  momentum  is  nearly  conserved.  Bird[9]  has  also  shown  that  the  non¬ 
conservation  of  the  angular  momentum  has  no  significant  effects  on  the  DSMC  solutions.  The  gas  flows  in  complex 
microfluidic  systems  are  expected  to  be  three-dimensional  and  unsteady.  As  an  eligible  direct  numerical  simulations 
tool  for  microflows,  it  is  then  important  for  DSMC  to  be  able  to  resolve  also  the  unsteady  flow  disturbances  with 
three-dimensional  structures.  DSMC  allows  for  a  direct  physical  simulation  of  gas  flows.  The  minimum  level  of 
physical  modeling  is  needed  at  the  level  of  molecular  models  that  have  been  used  in  kinetic  theory  of  gases. 
However,  a  large  number  of  particles  are  required  in  a  DSMC  simulation  to  reduce  the  inherent  statistical  scatter. 
The  computer  resource  needed  to  perform  three-dimensional  flow  simulations  can  be  significant,  which  appears  to 
be  a  major  limiting  factor  for  performing  many  necessary  simulations  to  capture  three-dimensional  physical  flow 
disturbances.  Many  efforts  have  been  made  to  take  advantage  of  the  advent  of  parallel  computing  to  reduce  the  time 
required  for  DSMC  simulations.  The  simulations  presented  here  have  been  performed  by  the  use  of  a  parallel  DSMC 
algorithm  [10]  and  a  parallel  computer  cluster,  CEPCOM,  in 
the  CFD  Lab.  at  Western  Michigan  University. 

In  this  paper,  results  of  three-dimensional  time- dependent 
microflow  simulations  using  parallel  DSMC  computation  will 
be  presented.  The  flow  geometry  considered  is  shown  in 
Figure  1 .  The  flow  develops  between  two  diffuse  walls  in  the 
y-direction,  the  top  wall  moves  at  a  speed  of  U.  The  lower 
wall  was  brought  to  a  stop  instantaneously  at  the  beginning  of 
the  calculations.  The  channel  height  h  is  nondimensionalized 
by  A  ,  the  mean  free  path  of  the  initial  gas.  Periodic  boundary 
conditions  are  applied  in  the  x-  and  z- directions.  An  artificial 
gravitational  forcing  (AF)  is  applied  along  the  y-direction  to 
the  molecules. 


AF  =  c(u-U)2 


(1) 


U  represents  the  fluid  velocity  component  in  the  x-direction.  c  denotes  a  constant  coefficient  that  is  used  to  adjust 
the  level  of  the  external  forcing  on  the  flow.  A  similar  flow  geometry  and  imposition  of  the  artificial  acceleration 
have  been  previously  developed  by  Malek  Mansour[ll]  in  an  analytical  study  of  the  onset  of  hydrodynamic 
instabilities  and  independently  by  Bird [12]  to  study  the  forced  micro  Couette  flows.  The  effect  of  the  artificial 
gravitational  force  is  that,  when  a  fluid  particle  is  perturbed  and  moves  in  the  y-direction,  an  effective  buoyancy 
force  favor  further  displacement  of  the  particle  in  the  y-direction.  The  scenario  is  similar  to  the  instability 
mechanism  associated  with  the  centrifugal  force  in  the  Taylor-Couette  and  Gortler  problems,  where  the  Taylor- 
Gortler  vortices  have  been  observed  in  experiments  and  computed  using  two-dimensional  DSMC.  By  using  an 
artificial  forcing  in  the  form  of  Equation  (1)  instead  of  a  physical  one,  it  is  possible  to  apply  a  wide  range  of  forcing 
level  to  the  flow  without  changing  the  flow  geometry.  Since  the  normalized  root  mean  square  of  the  disturbances  is 
inversely  proportional  to  the  square  root  of  the  sample  size, [13]  the  present  flow  model  allows  for  a  more  efficient 
use  of  the  computer  resources  in  capturing  unsteady,  three-dimensional  disturbances  by  artificially  energizing  the 
flow.  It  is  reasonable  to  expect  that  vortical  structures  or  flow  patterns  that  are  similar  to  the  Taylor- Gortler  vortices 
will  also  be  found  in  the  current  model  flow.  In  fact,  earlier  simulations[14,15]  of  the  same  flow  geometry  did 
capture  highly  organized  discernable  vortical  mean  flow  structures.  Fourier  analyses  of  the  vertical  velocity  signals 
show  the  existence  of  a  hierarchy  of  fluctuation  modes  that  consist  of  a  fundamental  and  its  harmonics.  Up  to  six 
harmonics  can  be  readily  identified.  In  the  following,  results  with  different  forcing  coefficient  will  be  presented  and 
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compared  for  two  different  channel  heights.  The  dynamics  and  the  flow  structural  pattern  of  the  sample  averaged 
flow  disturbances  generated  by  the  applied  forcing  will  be  examined. 


RESULTS  AND  DISCUSSION 


In  this  section,  results  of  the  simulations  will  be  presented.  The  box  height  is  320.  The  parallel  DSMC  code  used 
in  this  study  has  been  validated  by  comparing  with  the  serial  code  results  of  microchannel  flows  [2]  and  the  results 
were  reported.  [3]  For  the  present  flow  geometry,  a  comparison  has  also  been  made  with  the  serial  code  results  of 
Bird  for  a  2D  domain  that  slices  through  the  current  box  flow  along  a  z-plane.[14]  Flow  patterns  captured  by  the 
two  codes  were  basically  the  same. 

Two  types  of  sampling  have  been  used  in  the  study  to  obtain  macroscopic,  averaged  and  disturbance  components 
of  the  simulated  flows.  The  “unsteady”  properties,  say  for  the  vertical  velocity  component  v  is  obtained  by 

.  i 

v(x,y,z,t) - J  £  ^Mcvi(t  )dt  (2) 

*S  J  f~2ts 

where  cv  denotes  the  vertical  component  of  the  particle  velocity,  Mc  the  number  of  particle  in  the  cell  at  time  t' .  ts 
represents  the  short  time  period  for  the  unsteady  sampling.  It  was  set  at  twenty  times  the  mean  collision  time  of  the 
initial  gas.  The  long-time- averaged  or  averaged  component  of  the  flow  is  obtained  by 


—  f"?%  Old,' 

L~ tnJto  ^l=1  ’ 


(3) 


where  t0  denotes  the  beginning  time  of  the  sampling.  The  difference  between  the  unsteady  sample  average  and  the 
long  time  average  is  defined  as  the  perturbation  v' .  Or, 


v  =  v  -  v  (4) 

The  artificial  forcing  is  applied  in  the  three-dimensional  simulations  with  a  computational  box  of  size 
800x320x640.  The  box  dimensions  are  nondimensionalized  by  the  initial  mean  free  path  A  .  The  Knudsen  number, 
based  on  the  channel  height  h  is  1/320.  In  the  present  three-dimensional  simulation,  48  million  particles  were  used. 
Typically,  twenty  processors  were  used.  Depending  upon  the  level  of  forcing  and  the  box  sizes,  the  computing  time 
to  achieve,  for  instance,  stationary  results  are  regularly  in  the  order  of  two  to  three  weeks.  For  the  present 


simulations,  the  time  steps  are  20%  of  the  mean  molecular  collision  time  of  the  initial  gas  (Tc).  That  is,  At*  =  0.2 . 
In  a  previous  work, [15]  it  has  shown  that  unsteady  DSMC  simulations  are  independent  of  their  exact  values  as  long 
as  the  time  steps  are  small  compared  with  the  mean  molecular  collision  time.  It  is  valid  for  not  only  the  stationary 
state  of  a  flow  but  also  its  development  in  time.  The  speed  U  of  the  top  plate  is  set  equal  to  twice  the  most  probable 
molecular  speed  cmp  at  the  initial  gas  temperature  of  1,000  K.  That  is  U= 2.  The  Reynolds  number  is  1,156.  The 


forcing  level  is  changed  by  varying  the  value  of  c  (=0.001,  0.01).  The  maximum  forcing  that  can  occur  in  the  flow 
then  varies  by  a  factor  of  about  10.  The  results  for  the  same  flow  geometry  with  different  channel  heights,  plate 
speeds,  and  forcing  constants  have  been  shown  in  earlier  work. [15, 16]  For  comparisons  with  the  results  of  the 
current  simulations,  some  results  with  h= 60,  U=  1,  and  c=0.04  from  the  previous  simulations  are  shown  in  Figures  2 
and  3.  The  computational  box  size  is  800x60x320  and  16  million  particles  were  used.  The  Knudsen  number  is 
1/60.  Figure  2  shows  the  long-time-averaged  in-plane 
pathlines  at  three  axial  locations.  The  long  time  average 
begins  at  the  initiation  of  the  simulations.  The  mean  flow 
pattern  has  reached  a  stationary  state.  The  flow  is  periodic  in 
the  spanwise  direction  with  three  counter  rotating  vortical 
structures.  The  wavelength  is  320/3.  As  is  mentioned  earlier, 
the  vortical  structures  appear  as  a  result  of  the  externally 
applied  forcing  in  the  y-direction  that  mimics  the  centrifugal 
forcing  in  the  Taylor-Couette  problems.  Figure  3a  shows  the 
energy  spectra  of  the  instantaneous  (unsteady),  the  long-time- 
averaged,  and  the  disturbance  flow  fields.  The  velocity 
component  signals  are  taken  at  the  cell  that  is  located  at  the 
center  of  the  simulation  domain.  Note  that  a  disturbance  has 
been  defined  in  the  present  study  in  Equation  (4)  as  the 


FIGURE  2.  In-plane  pathlines  of  the  long-time- 
averaged  flow  for  h- 60,  U=  1,  and  c=0.04. 
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difference  between  the  unsteady  sample  average  and  the  long  time  average.  Figure  3b  shows  the  spectra  of  the  cross 
correlation  of  the  vertical  and  the  spanwise  velocity  components.  The  long-time-averaged  component  of  the 
instantaneous  velocity  field,  which  shows  highly  spatially  coherent  patterns  (Figure  2),  is  temporally  dominated  by 
the  zero  frequency  (=(m-l)/(20Tc)),  or  stationary  mode.  The  spectrum  distributions  at  other  frequencies  are  flat, 
indicating  that  the  signals  are  mainly  statistical  fluctuation  at  those  frequencies.  The  correlation  spectra  show  that 
the  mean  flow  is  highly  correlated  at  zero  frequency.  The  distributions  of  the  disturbance  correlation  spectra  are 
flat,  indicating  that  there  does  not  appear  to  be  very  strong  correlations  between  the  vertical  and  the  spanwise 
disturbance  velocity  components. 


- a - Instantaneous 

-f a —  Disturbance 


m 


- a - Instantaneous 

- ajjg|  -  Disturbance 


(a)  (b) 

FIGURE  3.  Fourier  spectra  for  h- 60,  U=  1,  and  c=0.04.  (a)  Energy;  (b)  Cross  correlation. 


With  the  same  flow  conditions  and  parameters  in  the  artificial  forcing,  more  work  is  done  to  the  particles  in  the 
larger  flow  domain  of  h= 320  than  those  of  h- 60  and  the  flows  are  more  energetic.  For  the  h- 320  cases  that  are 
shown  in  the  following,  the  speed  of  the  top  plate  U  is  set  equal  to  2  and  c=0.001  and  0.01,  respectively.  Figure  4a 

shows  the  in-plane  pathlines  of  the  mean  flow  at  three  x-locations  at  ^  ^  =83.08,  where  the  time  scale  h/U  is  about 
320  Tc.  The  mean  flow  pattern  has  reached  a  stationary  state.  The  periodic  flow  pattern  has  one  pair  of  vortical 
structures  in  the  mean  flow  in  the  spanwise  direction.  The  time  development  of  the  flow  shows  that  secondary  flows, 


FIGURE  4.  In-plane  pathlines.  c=0.001.  (a)  Long-Time- Averaged  flow;  (b)  Disturbances. 


which  did  not  exist  in  the  initial  conditions,  first  develop  from 
the  lower  region  of  the  flow,  where  the  applied  forcing  is  large. 
Flow  structures  with  a  large  range  of  length  scales  then  develop. 
At  a  later  time,  the  small-scale  flow  patterns  disappear.  In  a  time 
sequence  animation  of  the  computed  flow,  the  flow  structures 
were  seen  evolving  through  events  that  involved  the 
amalgamation,  generation,  and  destruction  of  vortical  structures, 
large  and  small.  Figure  4b  shows  the  in-plane  pathlines  for  the 
disturbance  flow  on  the  same  three  planes  in  the  x-direction  as 
those  in  Figure  4a.  Despite  the  statistical  fluctuations,  there  are 
identifiable  flow  patterns  that  are  vortical  and  similar  in 
appearance  with  those  shown  in  Figure  4a  for  the  mean  flow. 
Figure  5  shows  the  time  history  of  the  disturbance  velocity 
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FIGURE  5.  Evolutions  of  disturbance  velocity 
components  for  c=0.001. 
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component  signals  at  the  cell  that  is  located  at  the  center  of  the  simulation  domain  from  the  beginning  of  the 
simulation.  The  DSMC  disturbance  velocity  signals  contain  statistical  fluctuations  about  averages  that  appear  to  be 
stabilizing  toward  the  end  of  the  simulation.  Figure  6a  shows  the  energy  spectra  of  the  instantaneous  (unsteady),  the 
long-time-averaged,  and  the  disturbance  signals  (shown  in  Figure  5).  Figure  6b  shows  the  spectra  of  the  cross 
correlation  of  the  vertical  and  the  spanwise  velocity  components.  The  long-time-averaged  component  of  the 
instantaneous  velocity  field,  which  shows  highly  spatially  coherent  patterns  (Figure  4a),  is  temporally  dominated  by 
the  stationary  mode.  The  disturbance  spectra  contain  regions  of  low  to  medium  frequency  with  Fourier  amplitudes 
that  are  at  least  two  orders  of  magnitude  higher  than  the  statistical  scatter  apparent  at  the  high  frequency  end  of  the 
spectra.  The  long-time-averaged  and  the  disturbance  kinetic  energy  spectra  increase  nearly  linearly  (in  the  log  scale) 
toward  the  low  frequency  end.  The  correlation  spectra  at  the  low  frequency  are  somewhat  higher  than  the  other 
modes  that  appear  quite  random. 


(a)  (b) 

FIGURE  6.  Fourier  spectra  for  c=0.001.  (a)  Energy;  (b)  Cross  correlation. 


Simulations  with  c=0.01,  which  is  ten  times  of  that  for  the  previous  case,  have  also  been  performed.  The  larger 
forcing  further  energizes  the  flow.  Figure  7a  shows  the  in-plane  pathlines  of  the  flow  disturbances.  The  disturbance 
flow  pathlines  are  rather  smooth,  indicating  that  the  resolved  disturbances  are  not  dominated  by  statistical 
fluctuations.  The  flow  disturbances  are  vortical  and  are  highly  three-dimensional.  The  sizes  of  the  three- 
dimensional  vortical  structures  vary  from  as  large  as  the  channel  height  to  very  small.  As  the  simulation  proceeds, 
these  structures  go  through  mixed  series  of  growth,  amalgamation,  and  decay  processes.  Note  that  the  secondary 
flow  in  the  mean  (not  shown)  and  the  three-dimensional  disturbances  do  not  exist  in  the  initial  conditions  and  the 
flow  is  not  preferentially  perturbed  by  any  means  in  the  DSMC  simulations.  Figure  7b  shows  the  velocity 


FIGURE  7.  Flow  disturbances  for  c=0.01.  (a)  In-Plane  pathlines;  (b)  Velocity  component  signals. 


disturbance  signals  registered  at  the  center  cell  of  the  simulation  domain.  Compared  to  the  disturbance  signals 
shown  in  Figure  5  for  c=0.001,  the  disturbance  velocity  components  are  larger  in  magnitude.  The  frequent 
occurrence  of  these  large  and  small  amplitude  excursion  reflects  the  evolution  process  of  the  disturbance  flow 
structures  given  in  Figure  7a.  Figure  8  shows  the  kinetic  energy  and  the  correlation  spectra  from  the  velocity  signal 
recorded  at  the  center  cell  of  the  computational  domain.  For  the  mean  flow,  the  spectrum  for  the  zero-frequency 
mode  is  significantly  higher  than  all  the  other  modes.  The  disturbance  kinetic  energy  spectra  are  high  at  the  low  and 
mid  frequency  range  and  drop  off  by  more  than  one  order  of  magnitude  at  the  high  frequency  end.  Similar  drop-off 
in  the  spectra  from  the  mid  to  high  frequency  modes  by  about  one  order  of  magnitude  is  also  apparent  in  the  velocity 
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correlation  spectra  shown  in  Figure  8b.  These  spectra  show  that  the  energetic  modes  of  disturbances  in  the  low  to 
mid  frequency  correlate  well  and  form  the  three-dimensional  vortical  structures  shown  in  Figure  7a. 


CONCLUDING  REMARKS 

In  addition  to  the  mean  flow,  the  unsteady  DSMC  simulations  of  the  micro  Couette  flow,  with  sufficient  artificial 
forcing,  show  clearly  unsteady,  three-dimensional  flow  disturbances  that  are  coherent  and  well  correlated  both 
spatially  and  temporally.  Since  an  artificial  body  force  is  applied  to  energize  the  system,  the  captured  events  may 
not  be  realistic.  Nevertheless,  it  is  evident  from  the  present  results  that  the  DSMC  method  is  capable  of  resolving 
oscillatory  three-dimensional  disturbances  of  physical  significance.  In  a  recent  work,  we  examine  the  micro 
Rayleigh-Benard  problem  where  a  physical  buoyancy  force  naturally  arises  due  to  the  temperature  gradient  and  no 
artificial  forcing  is  used.  The  preliminary  results  of  the  unsteady,  three-dimensional  DSMC  simulations [17]  show 
Rayleigh-Benard  cells  in  the  mean  flow  and  disturbances  that  are  unsteady  and  three-dimensional.  These  findings 
suggest  that  unsteady  DSMC  can  be  used  for  the  direct  numerical  simulations  of  the  gas  flows  in  microfluidic 
devices. 
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